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In this paper, we initiate a new study of steady funnel-flow accretion onto strongly 
magnetized neutron stars, including a full treatment of shock generation. As a first step, 
we adopt a simplified model considering the flow within Newtonian theory and neglecting 
radiative pressure and cooling. The flow is taken to start from an accretion disc and then 
to follow magnetic field lines, forming a transonic funnel flow onto the magnetic poles. 
A standing shock occurs at a certain point in the flow and beyond this material accretes 
subsonically onto the star with high pressure and density. We calculate the location of 
the standing shock and all other features of the flow within the assumptions of our model. 
Applications to observed X-ray pulsars are discussed. 



§1. Introduction 



> 

in 

0^ ■ Following the discovery of periodically variable X-ray sources^ 1 ™ and subsequent 

theoretical explanations, systems consisting of strongly magnetized neutron stars 
accreting from their binary companions have become widely accepted as providing 
the model for periodic X-ray sourcesP 1 '^'^ The accreted matter coming via Roche 
lobe overflow from the companion forms an accretion disc around the neutron star 
(NS). If the NS has a strong magnetic field (^> 10 12 G), eventually at least some of the 
accreting matter cannot continue to accrete directly onto the NS surface in the disc 
plane but rather is constrained to follow the magnetic field lines, forming a funnel 
flow onto the magnetic poles (e.g. Ref. U])). The polar caps reach high temperatures, 
emitting high-energy radiation in the form of X-rays. If the magnetic axis and the 
rotational axis are misaligned, this high energy radiation from the magnetic poles 
may be observed as periodic pulsations PJSKELE} 

Polar accretion models assuming strong magnetic fields have been studied by 
many authorsP'®'^' 1 ^^'^'^'^' 1 ^ In these studies, it has been shown that ac- 
cretion onto the polar regions of strongly magnetized NSs and shock generation near 
to the NS surface can explain the high energy emission seen in X-ray pulsars. Al- 
most all previous studies, however, have been limited to cylindrical or conical funnel 
flow aligned with the magnetic axis, whereas the actual accretion flow must have 
a curved geometry along the magnetic field lines going from the disc to the stellar 
pole. Including the effects of the curved geometry may well be very important for 
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understanding the complex behaviour of the light curves of X-ray pulsars P' 1 ^ We 
here solve the flow equations in a consistent way, starting from where the funnel-flow 
material leaves the plane of the accretion disc and following it until it impacts on 
the magnetic poles. 

There have been several previous studies of this type which have included the 
effects of the curved flow geometry.^ JIDiEDfilD All of these solve the Bernoulli equa- 
tion, obtaining stationary transonic flow solutions similar to Bondi accretion!^ In 
the present study, we mainly follow the approach of Koldoba et al. (2002, hereafter, 
KLUR)P3 (We note that Koldoba and collaborators have subsequently carried out 
very interesting related calculations using an MHD computer code® but these were 
not addressing the question of shock location which is our main interest here.) As- 
suming a dipole magnetic field whose axis coincides with the stellar rotation axis, 
KLUR constructed the equations for transonic flow and calculated accretion flows 
which accelerate from the subsonic regime into a supersonic regime. In their study, 
however, although they found the location of the transonic critical point, they did 
not calculate the position of the shock front. In the present work, we now extend 
this, following the same basic picture but calculating also the shock location and 
strength and the effects of shock heating. This represents the first step in a larger 
project to study this mechanism in full detail. 

We adopt here a simplified model, considering a one-dimensional steady funnel 
flow along the field lines of a strictly dipole magnetic field aligned with the rotation 
axis, treating it within Newtonian theory and neglecting radiative pressure and cool- 
ing. Since the field is taken to be a purely dipole one (in the region of the funnel 
flow) we do not have any Poynting flux term in our calculation. We require the field 
to be strong enough to keep the funnel "rigid" and to keep the fluid flowing along the 
field lines: under the circumstances envisaged, this condition should be well-satisfied 
everywhere apart from close to where the flow leaves the disc and where it impacts 
on the NS surface. 

We recognize that the simplifications which we are making are very significant 
ones but they give us a problem which can then be solved completely in a straight- 
forward way, providing a useful basis for future work. Within the picture which we 
are using here, we can calculate the properties of the flow all the way from where 
it leaves the plane of the disc up to the stellar surface, allowing us to see several 
aspects which have not been observed in previous studies. In particular, we have 
shown that the shock location is sensitively dependent on the ratio of the pressures 
at the points where the flow leaves the disc and where it impacts on the NS surface. 
This is a crucial property of the accreting systems. Determining the shock location 
and the global characteristics of the flow is important for linking the properties of 
observed X-ray pulsars to the physical conditions in the NS - disc systems. 

We want to stress that our approach here does not attempt to include any 
detailed treatment of the physics near to the neutron star surface or near to where 
the funnel flow moves out of the plane of the disc. This has been treated elsewhere 
by a number of authors (as described in a later section). We focus instead on how 
the global solution fits together, something which we consider as certainly forming 
a necessary framework within which to fit these other more detailed discussions of 
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particular parts of the flow. 

The plan of this paper is as follows: in Section [2] we describe how we set up the 
problem. Section [3] then introduces the formulation for our treatment of transonic 
flow along the dipole magnetic field lines and Section U] presents the equations for 
treating the standing shock. Section [5] describes our method of solution for flows 
including shocks and presents our results. In Section El several related issues are 
discussed. Finally, Section [7] contains a summary of the paper together with some 
further discussion of the results and their astronomical implications. Throughout 
the paper, we use standard spherical coordinates (r, 9, tp) as well as cylindrical coor- 
dinates (R,(p,z), which are convenient for some of the descriptions. The subscripts 
"NS" and "d" denote, respectively, that the values of the physical parameters con- 
cerned are measured at the NS surface or at the point where the funnel flow leaves 
the plane of the disc. 



§2. Problem set-up 



We are considering here a strongly magnetized neutron star with a dipole field 
surrounded by an accretion disc. At some point, moving in through the disc, the 
accretion flow becomes significantly altered due to its interaction with the magnetic 
field which, in turn, is distorted away from its purely poloidal form within the disc 
because of interaction with the disc material (see Kluzniak Sz Rappaport 201 for a 
recent treatment). Eventually, at least some part of the accreting material starts 
to go along the dipole field lines onto the magnetic poles of the star. As with the 
previous studies, we consider here a rigidly rotating dipole field everywhere within 
the region of the funnel flow, co-rotating with the neutron starP^n^'D^'E^ For a 
sufficiently strong magnetic field, each element of the accreting matter will fall onto 
the star along a single field line passing through its initial position on the disc. The 
geometry of the dipole field can be described by 

r = i? d sin 2 6>, (1) 

where, is the distance from the centre of the star to the point where the funnel 
flow leaves the plane of the disc. 

The dipole field strength can be written as 



^ = 5r(4-^H, (2) 



r 6 V Rd 



where B p is the poloidal magnetic field and /x is the dipole moment. This equation 
gives the magnetic field strength as a function of radial coordinate r along a field line 
which intersects the disc at r = Rd- The location of the point at which the funnel 
flow leaves the disc is determined by the competition between pressure, gravity and 
magnetic forces and is hard to estimate precisely, particularly in view of instabilities 
which play a role in this. We therefore simply choose suitable values for R^ and 
carry out our calculations with these (see Section 5). 

Following KLUR, we note that transonic solutions are only possible for R^ close 
to the co-rotation radius R CO rot (i-e., where the disc is co-rotating with the neutron 
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star) and so we make the approximation of assuming that the stellar angular velocity 
Q coincides with the Keplerian rotation rate at r = R^; i.e., 

where G is the gravitational constant, and m is the mass of the star. 

Under these conditions, the transonic accretion flow can be solved for using the 
method given by KLUR (see Section [3|) . This flow will have a shock, at a certain 
point in the supersonic region, which we here assume to be adiabatic and planar. 
Next, we introduce the further assumption that there are no radiative energy losses 
from the inflowing material; the Bernoulli equation can then be applied everywhere 
along the flow-lines. 

Following these preparations, we can then solve for the whole stream configura- 
tion, including an adiabatic shock, if we impose suitable boundary conditions at the 
point where the funnel flow leaves the plane of the disc and at the stellar surface. 
The main goal of the calculation is to estimate the position of the standing shock 
and the distribution of physical quantities in the flow. 



§3. Formalism for flow without a shock 



Following KLUR, we solve the energy conservation equation (Bernoulli equation) 
with the condition that each fluid element moves along a magnetic field line going 
from r = i?a on the disc to the NS surface. For the present situation, we write the 
Bernoulli equation (in the reference frame co-rotating with the field line) as 

E=-v 2 + fiVsn 2 fl, 4) 

2 7- 1 r 2 K J 

where E is the specific energy (assumed to be constant), v is the poloidal velocity 
along the field line, a is the sound speed and 7 is the adiabatic index. The last 
term in this plays the role of a potential corresponding to the centrifugal force. 
Since we are considering a rigidly rotating dipole field throughout the region of the 
funnel flow there is, by assumption, no toroidal field component present there. This 
assumption (which is a rather reasonable one for most of the funnel flow) allows us 
to omit the Poynting flux term which would otherwise appear in Eq. ([4]). The role 
of the magnetic field then consists only in fixing the configuration of the flux tube, 
in keeping the fluid elements moving on it and in determining the variation with 
position of pv via the MHD mass conservation equation 

Anp^- = K, (5) 

(where K is a constant).^ 1 By using Eqs. ([2]) and ([5]), the kinetic energy of the fluid 
related to motion along the field line can be rewritten as v 2 jl = K 2 B 2 /32ir 2 p 2 and 
inserting this into the Bernoulli equation gives 

/ ^ K V ( A 3r\ sp<~ x 
32iT z p z r \ Rd J 7 — 1 
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Here we have used the specific entropy, s = a 2 /p 7_1 instead of the sound speed. With 
the assumption that the specific energy E and the specific entropy s are conserved 
along the flow-line, this equation gives a relationship between p and r and so we 
can calculate the density distribution p{r) along the flow-line if we specify boundary 
conditions at the point of departure from the disc: E = E(R^) = const, and s = 
s(Rd) = const. 

For convenience of calculation, we rewrite Eq. (|6|) in the dimensionless form: 
~ , , 4-3f sp- 1 (\ r 3 \ 

E = e(r, P) = ™r + T ~ ~ + ■ ( 7 ) 



2p 2 f 6 7-1 \f 2 

Here, we have used following transformation of variables, 

tj ~ Kp . 

r = Rdr > P= 4^Rj P > 

s = ^d(^^) 7 ^~ E = tfRlE. (8) 

(From here on, since we use only dimensionless quantities, we will omit the tildes.) 
By solving Eq. (|7|), we can obtain the values of the physical quantities everywhere 
along a flow-line if there is no shock. Note that when working in terms of the 
dimensionless quantities, the field strength does not itself enter the calculation, only 
its variation with position. 

We are interested in transonic solutions with the initially subsonic inflow becom- 
ing supersonic at a critical point. In order that the flow should be non-singular at the 
critical point (as required for a physically meaningful steady-flow solution) regularity 
conditions must be satisfied there, as given by Eqs. (16) and (17) of KLUR. 

§4. The treatment of shocks 

We will denote the location of the critical point by r c ; a shock may then appear 
beyond this in the supersonic part of the flow with the location depending on the 
boundary conditions. We next summarize our methodology for treating the shock 
within our simplifying assumptions of considering a stationary planar shock which is 
adiabatic (i.e. no energy is lost from the flow)P5J From the conservation of the mass, 
momentum and energy fluxes across the shock, one obtains the Rankine-Hugoniot 
relations linking the values of quantities on either side of it. The ones of relevance 
for us here are: 



(t + i)m? 

(7 _l )M 2 + 2 ' 



P2 = - 1W2 l , n Pl (9) 



2 7 M? - (7 - 1) 

7 + 1 
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where M = v/a is the Mach number and the subscripts 1 and 2 denote the values of 
flow quantities immediately before and after the shock. The shock raises the specific 
entropy of the fluid, with the change being given by 



In the next section, we use these shock relations together with the flow equation ((Tj). 



5.1. Flow beyond the shock 

The behaviour of the fluid quantities along the flowdine, for given boundary 
conditions at r = R^, is completely described by the Bernoulli equation (|7|) except 
that, when passing across the shock front, the value of the specific entropy needs to 
be changed to the new higher value, i.e. beyond the shock, the Bernoulli equation 
becomes: 



where S2 is the specific entropy of the flow beyond the shock, given by Eq. (jlip . 

Our procedure for calculating the entire flow solution with the shock, starting 
from conditions at r = R,x on the disc and ending at the surface of the NS, is as 
follows: (i) Choose a position for the shock r = r s , with tns < r s < r c . (ii) Solve 
for the flow in the region upstream of the shock, r s < r < using Eq. (|T|). This 
gives the value of pi = p(r s ), just before the shock, from which can be calculated the 
corresponding values of p\ and M\. (hi) Using the Rankine-Hugoniot relations, we 
then calculate the values of density and pressure beyond the shock [pi and P2) and 
hence the new value of the specific entropy S2- (iv) Using this value of S2, we then 
calculate the flow beyond the shock using Eq. (|12|) . In this procedure we obtain the 
value of p at the stellar surface and hence also the values of the other fluid parameters 
there, (v) In practice, we would prefer to specify the value of the pressure at the 
stellar surface rather than the location of the shock and so we iterate the solution 
until the shock location corresponding to the desired surface pressure is found. 

In Fig. (TJ we show our solutions for density and velocity in the special case 
of transonic flow with no shock. We have made calculations for 7 = 4/3 with 
parameter values of s = 0.003 and E = —1.448. Note that 7 = 4/3 is smaller 
than the maximum value of 7 for which transonic flow can occur in Bondi accretion 
(7 = 5/3). This figure also shows the increased entropy S2 which would appear 
beyond shocks occurring at different locations r s along the flowdines. The value 
of S2 is one of the most important quantities here, since it controls the whole flow 
solution in the downstream region, r^s < T < r s . 

Figs. [2] and [3] show the variation of density and pressure along the flowdines; the 
lowest solid curve corresponds to transonic flow without a shock, as shown in Fig. [H 
and the upper curves are for flows having an adiabatic shock at different locations. 
When there is a shock, density and temperature have discontinuous jumps across it, 




(11) 



§5. Accretion flow with a shock 




(12) 
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Fig. 1. Density and poloidal velocity distributions along flow- lines with no shock for different values 
of 7. We show the case with 7 = 4/3. Behaviour of the the enhanced entropy which would 
appear beyond shocks occurring at different locations r s is also shown. 

as given by Eqs. ([9]) and (jlOp . The envelope of the density values on the downstream 
side of the shock, p2, is shown by the dotted curve in Fig. [21 and the envelope for 
P2 is shown in the same way in Fig. The behaviour of the physical quantities for 
r NS < r < r s is given by Eq. (|12p . Note that the density can increase across the 
shock by at most a factor of (7 + 1)/ (7 — 1), within the non-relativistic theory being 
used here, whereas increases of several orders of magnitude may occur by means of 
progressive compression in the subsonic flow beyond the shock. 

From Figs. [2] and O we see that the resulting density and pressure at the stellar 
surface depend sensitively on the location of the shock (as an example, tns = 0.005i?d 
is shown as the vertical dotted line in the figures). As mentioned previously, if the 
pressure (or density) at the stellar surface is given as a boundary condition, the 
shock location can then be determined and the whole solution completed. Note that 
in Figs. [2] and El the curves for the higher values of r s rise quickly just beyond the 
shock and then flatten to an almost constant slope, whereas those for smaller r s 
have an almost constant slope throughout the region beyond the shock. The reason 
for this difference can be understood from Eq. f)12|) : in calculating the behaviour of 
the density and pressure just beyond the shock, the gravity term oc r _1 completely 
dominates over the centrifugal term oc r 3 when the shock occurs near to the stellar 
surface, whereas the centrifugal term makes a significant contribution with respect 
to the gravity term further from the star. This difference is reflected in the two- 
component gradients of the curves in Fig. U introduced below. 

5.2. Pressures and densities achievable at the NS surface 

Before going into details about the shock location, we here discuss the range of 
possible densities in the flow at the point where it impacts on the NS surface. The 
maximum achievable density corresponds to the shock occurring just after the flow 
passes the critical point. In this case, the fluid reaches the NS surface with almost 
zero poloidal velocity and so the kinetic energy term in Eq. (|12p can be neglected 
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Fig. 2. Density distributions for flows including shocks, taking 7 = 4/3 and s = 0.003. The upper 
solid curves are for flows with shock discontinuities at r = 0.85,0.5,0.1 and 0.025 (with the 
dashed lines marking the shock discontinuities); the dotted curve is the envelope of the density 
values on the downstream side of the shock; the lowest thick solid curve is for a flow with 
no shock. The vertical thin-dotted line denotes the position of the NS surface in the case of 
rNs/Rd = 5.0 x 10" 3 . 




Fig. 3. Pressure distributions corresponding to the density distributions in Fig. [2] 

theraEl The rotational kinetic energy term at the NS surface can also be neglected 
since the field lines are taken to be co-rotating with the neutron star which is itself 



*' This can be easily understood if we consider the analogy with Bondi accretion. Bondi flow 
has two trans-sonic sequences, one representing flow which is being accelerated from sub-sonic to 
super-sonic and the other representing flow being decelerated from super-sonic to sub-sonic. An 
accretion-flow solution starts off along the first (accelerated) sequence but, after the critical point, 
a shock may occur in which case the solution jumps down to the second (decelerated) sequence, 
where it is then sub-sonic. Following this, the flow is progressively decelerated. If the shock occurs 
very close to the critical point, there is a relatively long deceleration time before the flow reaches 
the stellar surface and so the velocity there will be quite low. 
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slowly-rotating for all of the cases of interest here. We then obtain 



7-1 




) 



1 1/(7-1) 



/>NS,max 



(13) 



The minimum achievable density corresponds to the shock occurring just above the 
stellar surface. In this case, the kinetic energy dominates over the thermal energy 
just before the shock and we obtain 



The first term on the right hand side is the factor appearing in Eq. ([9]). 

For given outer boundary conditions (E, si, R d ), the density in the flow at the 
NS surface must lie within the range between these maximum and minimum values 
for any flow including a shock. The range of possible pressures in the flow at the NS 
surface can be obtained in the same way. 

5.3. Location of the shock 

Figs. [2] and [3] show values for the normalized density and pressure along the 
flow- lines, plotted against the normalized radial coordinate for various shock loca- 
tions. We take a canonical value of 10 km for the neutron star radius but this then 
corresponds to different values of r depending on the value taken for R d , the lo- 
cation at which the funnel flow leaves the disc. (For our calculations, we use the 
following representative values for the other neutron star parameters: m = 1.4M Q , 
li = 10 30 G cm 3 giving -B P) ns ~ 2.0 x 10 12 G - see, for example, Shapiro and Teukol- 



At the outer and inner boundaries of the funnel flow we impose pressure-matching 
boundary conditions. When we know the starting point of the flow (Rd) and the 
ratio of the pressures at the inner and outer boundaries {pNs/Pd), we can obtain the 
appropriate shock position. Fig. U]shows the shock position r s , normalized by the NS 
radius tns, plotted against the pressure ratio PNs/.Pd for a range of selected values 
of R d (r NS /R d = 0.1,0.05,0.02,0.01,0.005,0.002 and 0.001, as marked next to the 
related curves), taking account of the fact that the actual values of R d are rather 
uncertain. The shock positions (in a log scale) are distributed within a band roughly 
proportional to log(pNs/Pd): as the pressure ratio becomes larger, the shock moves 
further away from the NS surface. The shock positions obtained using the various 
specific boundary conditions discussed below are indicated by the filled circles in this 
figure. We now turn to discussion of these conditions. 

5.3.1. Outer boundary condition 

At the outer boundary r = R d , we set the pressure in the funnel flow equal to the 
gas pressure at that point in the disc. We discuss here a number of possible models 
for this although all of them should be considered as being very approximate. 

First, we consider the picture where the outer boundary of the funnel flow is 
matched onto a Shakura-Sunyaev discP^ with the funnel flow being taken to start 



( 7 + 1)M? >NS r R^R d - 3r NS ) 
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Fig. 4. Shock locations as functions of the pressure ratio. The filled circles correspond to the cases 
listed in Table|T]as indicated by the adjacent labels. The trend for the shock moving further out 
as the pressure ratio increases was to be expected since the stopping power of the back pressure 
becomes greater as the pressure ratio increases, causing it to be felt further out. The reason for 
the change in slope seen within each individual curve has been discussed in Section \5. II 



from the point where the magnetic pressure of the dipole field balances the gas 
pressure.® This gives 

Pd = Bp [ Rd)2 = 4.65 x 10 15 a- 9 / 10 

where a is the alpha viscosity parameter (here we restrict attention to the situation 
for a = 0.1). Using B P (R&) 2 /Si: = fi 2 and fixing m cr it, this gives the value of 
R& (™crit is the critical accretion rate which gives the Eddington luminosity). We 
consider models with two values of m/m cr it. "Model A" with m/m cr it = 2/3 and 
"Model B" with m/m crit = 1. 

Our third model ("Model C") is based on the disc model of Kluzniak & Rap- 
paportpSJ w ho derived a modified form of the Shakura-Sunyaev solution for discs 
around magnetized stars, taking account of the way in which the magnetic field 
modifies the flow in the inner regions of the disc and is, in turn, modified by its in- 
teraction with the flow. They used two forms of ansatz for the toroidal component of 
the field produced in the disc (we follow their first prescription here) , and calculated 
the effect on the flow of the resulting magnetic torque which progressively substitutes 
the standard viscous torque as one moves in through the disc. This gives a set of 
formulae for the disc profile and the variations of the fluid parameters through the 
disc which represent generalisations of the corresponding Shakura-Sunyaev formulae. 
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The standard Keplerian disc then ends at the torque balance point r = i?torq (where 
the magnetic torque has grown to balance the accretion torque), i.e. where 

-^torq-^Pitorq-S^torq = Tfl 

For the parameter values mentioned above, this comes very close to the co-rotation 
radius R CO rot if the stellar rotation period is at the lower end of the range for X-ray 
pulsars. We again take the funnel flow to leave the plane of the disc at the point 
of balance between the gas pressure and the magnetic pressure of the dipole field; 
this balance occurs very close to the torque balance point. Our Model C has a 
Kluzniak-Rappaport disc with the same overall parameters as Model A above and 
has an NS rotation period of 1 second for which i?torq is indeed very close to R CO rot 
(and i?d ~ -Rcorot as required for our overall picture). 

Finally, we consider the situation when the outer accretion flow consists of matter 
coming essentially in free-fall from a companion star and the funnel flow starts at 
the Alfven radius r& given by 

|f = \ P {TA)v 2 {r A ). (17) 

We then take = ta and we refer to this case as our "Model D" . 
5.3.2. Inner boundary condition 

At the inner boundary r = tns; the fluid flow changes from being essentially 
radial to being transverse and the situation there is even less clear than at the outer 
boundary since we must consider a matching with the transverse pressure (including 
magnetic pressure) at a suitable height in the NS atmosphere. The flow pressure at 
the NS surface might then be as high as that given by balance with the magnetic 
pressure pns = -^pNs/^ 71 " ~ 10 23 dyn cm~ 2 . 

Sophisticated NS atmospheric models give various different views about the NS 
surface pressure. For instance, Van Ripei® found that the density of the outermost 
layer of the atmosphere of a magnetized NS should be in the range ~ 0.5 — lOg cm~ 3 , 
for m pa l.OM , r^g pa 10 6 cm, and -Bns ~ 10 12 G. However, the corresponding 
pressure is strongly temperature dependent and so even if we fix on p pa lg cnr 3 , 
the pressure can take a wide range of values going from ~ 10 16 dyn cm -2 up to 
~ 10 22 dyn cm~ 2 for the temperature ranging from T ~ 10 7 K to ~ 10 9 K. The 
temperature at r = tns is almost independent of the shock position (as discussed in 
Section 5.4), giving T NS = 2.4 x 10 9 K for a = 0.1. According to Van Riperj^ for 
Tns ~ 10 9 K, the pressure of the NS surface layer is roughly ~ 10 21 — 10 22 dyn cm -2 . 

On the other hand, Brown and Bildsten 291 suggested that the accreted matter 
cannot spread out until it reaches a greater depth, since the magnetic pressure is 
always much larger than the matter pressure at the polar cap. In this case, our 
pressure boundary condition should be modified to ~ 10 25 [dyn cm -2 ], and the height 
of the shock above the neutron-star surface would become much larger. 

Another example of an atmosphere model for an accreting NS is the plane- 
parallel thin-slab model of Harding et al! 2 -^ According to this, the density in the 
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atmosphere changes by about five orders of magnitude within a thin layer of width 
<^ 200 cm with very little accompanying change in temperature. 

Detailed consideration of the physics of the near-surface region is beyond the 
scope of the present paper which focuses instead on global properties of the flow. 
Because of the uncertainties involved, we consider here a range of values for p^s, 
bearing in mind the atmospheric pressures quoted in above references, and we will 
return to give a more complete treatment of the near-surface region in future work. 

From Fig. HJ it can be seen that the shock location comes closer to the NS 
surface when either (1) the NS surface pressure is low, or (2) the pressure is high at 
the point where the funnel flow leaves the disc, i.e. when the ratio PNs/Pd is small, 
the shock will occur near to the NS surface. In the opposite limit, the shock will 
occur far above the NS surface, similar to the situation considered by Basko and 
SunyaevP As discussed in Sec. 16.21 below, the steep pressure gradient occurring in 
an optically thick radiative accretion column would effectively lead to a rather low 
NS surface pressure, bringing the shock nearer to the NS surface. Results for the 
various specific disc models, and for a relevant range of values of pns> are listed in 
Table U and depicted in Fig. SJ Cases AO, BO, CO and DO are the ones with the 
optically thick radiative columns. 

Table I. Shock positions for selected surface conditions: 
Shakura-Sunyaev disc with moderate accretion (A), 
Shakura-Sunyaev disc with critical accretion (B), 
Khizniak-Rappaport disc with moderate accretion 
(C) and the Free-fall model (D). The cases with opti- 
cally thick radiative accretion columns are indicated 
by a *, and for those models, the inner boundary 
condition for the funnel flow is imposed at the top of 
the column (see Sec. I6.2p . 



Disc Model 


NS Surface Condition 


Shock Position 


Label 




p N s[dyn/cm 2 ] 


r B [cm] 




Shakura-Sunyaev disc with 


1.0 x 10 i8 (*) 


1.9 x 10 b 


AO 


m/rhcrit = 2/3 


1.0 x 10 19 


3.9 x 10 6 


Al 


R d = 4.4 x 10 7 [cm] 


1.0 x 10 21 


25.0 x 10 6 


A2 


p d = 5.6 x 10 12 [dyn/cm 2 ] 


1.0 x 10 23 


32.9 x 10 6 


A3 


Shakura-Sunyaev disc with 


1.0 x 10 i8 (*) 


2.0 x 10 fa 


B0 


m/?Ti cr it = 1 


1.0 x 10 19 


2.6 x 10 6 


Bl 


R d = 3.7 x 10 7 [cm] 


1.0 x 10 21 


15.6 x 10 6 


B2 


p d = 1.5 x 10 13 [dyn/cm 2 ] 


1.0 x 10 23 


25.9 x 10 6 


B3 


Kluzniak-Rappaport disc 


1.0 x 10 i8 (*) 


18.1 x 10 B 


CO 


rh/rhcrit = 2/3 


1.0 x 10 19 


50.9 x 10 6 


CI 


Rd = 1.6 x 10 8 [cm] 


1.0 x 10 21 


108.4 x 10 6 


C2 


p d = 2.0 x 10 9 [dyn/cm 2 ] 


1.0 x 10 23 


136.3 x 10 6 


C3 


Free-fall model 


1.0 x 10 18 (*) 


22.4 x 10 e 


DO 


Rd = 3.1 X 10 8 [cm] 








p d = 4.6 x 10 7 [dyn/cm 2 ] 
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§6. Discussion of related issues 

6.1. Departure of the funnel flow from the plane of the disc 

In discussing the disc models above, we estimated the location of the point where 
the funnel flow leaves the plane of the disc by assuming that the magnetic pressure of 
the dipole field of the NS balances the pressure of the accreting gas there, following 
le & ReesP This seems a reasonable approximation but there are a number 
of issues that make the situation less clear. First, we should stress that there is a 
distinction to be made between the place where material starts to leave the plane 
of the disc and the place where one should define the "inner edge" of the disc to 
be, although they are often treated as being the same. In an accreting system, the 
"inner edge" is not a completely clear concept although one can take it to be the 
place where the disc ceases to be nearly Keplerian in the case of a geometrically thin 
disc. Also, while one often takes the disc to be "thin", nevertheless it does have a 
vertical structure and conditions are certainly not constant through the height of it. 
It is likely that there will be a separation between some material which goes into the 
funnel flow while other material continues in or near the disc plane. Another point 
is that the pressure-balance argument is based on supposing that the poloidal field 
of the neutron star is not affected by its interaction with the disc material whereas, 
in fact, there will be at least some region where the material flow distorts the field 
lines producing a toroidal component within the disc and giving a torque on the 
accreting material, as mentioned earlier. Finally, there is the issue of a variety of 
instability mechanisms which can influence the situation. It is therefore clear that 
the question of how and where material leaves the disc plane to go into the funnel 
flow is a complicated one. 

Following the pioneering study by Ghosh and Lamb, 111 many authors have 
worked to improve understanding of the magneto-hydrodynamic properties of ac- 
cretion discs but this still remains a difficult subject. Frequently, the position of 
the "inner edge" of the disc is taken to come at the point of balance between the 
magnetic torque and the accretion torque (see the discussion in Section 15.3. ip with 
the flow remaining essentially Keplerian down to that point. This is the case for the 
Kluzniak-Rappaport model which we have useo® and there one finds that the pres- 
sure balance point is extremely close to the torque balance point for parameter values 
of interest here. However, this depends on the ansatz made there for describing the 
toroidal field and, while that seems a rather reasonable one, this is something 
which should, in the end, come from a calculation including all of the physical fea- 
tures concerned. In basic MHD treatments of discs, the effective inner radius is 
almost coincident with the co-rotation radius, -Rcorotp^'^ but some more detailed 
studies suggest that it may be significantly inward of this, at -Rtorq ~ 0.1i? coro tP^ 
Alternatively, its location may be related to the Alfven radius, as in Model D in our 
treatment!^ 

Even if the gas pressure does balance the magnetic pressure at the starting 
point of the flow, the magnetic pressure would then grow very rapidly along the 
flowlines (as ~ r -6 ) and would soon become dominant so as to establish the funnel- 
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like geometry of the flow. This would then be expected to persist until the gas 
pressure could catch up with the magnetic pressure again after the flow penetrates 
into the NS surface layers. 

6.2. Radiative effects 

In the above treatment, we have not included the effects of radiation emitted 
by the fluid comprising the funnel flow but, in realistic situations, radiation pressure 
and radiative energy losses will be important near to the stellar surface. We will 
include these effects consistently in our future work but here we make just a rough 
estimate of how radiation pressure would change the present results. For doing 
this, we join our flow solution onto the solution for an accretion column calculated 
by InoueJ-^ which does include radiation pressure. Inoue (1975) solved the fluid 
equations, including photon diffusion, for a short accretion column just above the 
NS surface and found that the density distribution can be written as 



where, r^is = K^m/Sirc with being the Thomson scattering opacity. This den- 
sity distribution (and the corresponding pressure distribution) has quite a strong 
dependence on r and the accretion column becomes optically thin rather close to the 
NS surface. In this solution, the pressure of the column decreases by roughly five 
orders of magnitude within the optically thick regime. 

According to Inoue, the height of the optically thick radiative zone will be ~ 
6 x 10 cm for to = 10 17 g s" 1 . 101 When we join our solution at the top of the column 
with the density and pressure obtained from Eq. (|18p . the resulting shock position is 
closer to the stellar surface. Taking the pressure at the top of the radiative column 
(r = 1) to be ~ 10 18 dyn/cm 2 and the outer boundary condition to be that for 
Shakura-Sunyaev disc model with to = |m cr it — 10 17 g cm -1 (Model A), the shock 
position is then just above the top of the radiative column. This means that the flow 
becomes optically thick just after the shock and then emits high energy photons as 
black body radiation. If, on the other hand, we take the same pressure at the top of 
the radiative column but use the free-fall condition (Model D) for the inner boundary 
condition, the shock position is then still far from the NS surface and from the top of 
the radiative column (r s = 22.4 x 10 6 cm). Hence, in this case the flow becomes dense 
and achieves high temperatures beyond the shock, but it is still optically thin there 
and so emits energy via free-free radiation or cyclotron processes before forming the 
optically-thick radiative column. This optically-thin, but actively radiating region 
corresponds to a "transition region".^® 

6.3. Temperature distribution 

In this subsection, we present some results for the temperature distribution be- 
yond the shock, as estimated using T = S7 _1 p 7 . In Fig. [5j the temperature 
distribution for a transonic flow with no shock is shown by the lower solid curve 
and the corresponding temperature distributions for shocked flows with the selected 
values of r s are also shown (cf Figs. [2]and[3|). It is interesting that the tempera- 
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ture curves beyond the shock are independent of the shock position. This can be 
understood from the solution of the Bernoulli equation: since the kinetic energy is 
negligible beyond the shock, we have 



Toe 



7 



1 



7 



31 



I r 
E + - + — 
r 2 



(19) 



Since 7 and E are constants, the temperature behaves as a function only of r and 
does not depend on the shock position. 




0.01 0.1 1 



r 

Fig. 5. Temperature distribution along the flow-lines, taking 7 = 4/3. 



We note that the shock is not in the highest-temperature part of the flow (al- 
though the temperatures beyond the shock are considerably higher than those before 
it) and that there are parts of the shocked flow at relatively lower temperatures when 
r s becomes larger. Radiation from that part of the flow may appear in the low en- 
ergy components of the observed spectra of X-ray pulsars. Also, if the shock heating 
occurs away from the stellar rotational axis, relatively hot fluid will exist in a broad 
region covering the stellar pole. The distribution of hot matter may change the pulse 
shapes: if the accretion flow emits radiation from a broad region, the geometry of the 
emitting flow will be changed and this may change the shape of the pulses (see, for 
example, Ref. 16)). As noted above, the estimates for the shock location given here 
are not precise but a comparison between these shock locations and the observed 
properties of X-ray pulsars would be interesting. 

Fig. [5] shows temperatures measured in units of and for converting these into 
actual temperature estimates, one needs to use suitable models for the disc and for 
the departure of the funnel flow from it. We feel that it is necessary to be cautious 
about interpreting these results directly in terms of actual temperatures in the flow 
and at the neutron star surface, both because of the diversity of the models and also 
because of the simplified treatment in this sub-section. However, we note that surface 
temperatures derived from this, for the models which we have been considering, lie 
in the range from ~ 10 9 K (standard disc case; this value is consistent with values 
often quoted: see, for example, Ref. |29|) ) to ~ 10 11 K (free-fall case; this value is too 
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high for actual X-ray pulsars and would indicate that energy losses from the funnel 
flow needed to be taken into account). 

§7. Summary and further discussion 

In this study, we have considered funnel-flow accretion onto a magnetized neu- 
tron star (NS), including taking into account the possible occurrence of a standing 
shock. We have assumed that the funnel flow goes strictly along the magnetic field 
lines and that the field is a dipole field, aligned with the stellar rotation axis, through- 
out the region of our calculation. The funnel flow leaves from the plane of the disc 
and accelerates smoothly into the supersonic regime, passing through a critical point. 
Then, at some point after this, a standing shock occurs beyond which the flow then 
accretes sub-sonically onto the neutron star. We have calculated the location of 
the standing shock, neglecting energy losses, and have obtained a full solution for 
given boundary conditions at the stellar surface and at the point where the funnel 
flow leaves the disc. The results suggest that the shock position depends strongly 
on the boundary conditions. As the ratio between the pressures at the inner and 
outer boundaries becomes larger, the shock location goes further from the NS sur- 
face. If the NS boundary pressure is high and free-fall type accretion occurs, the 
shock location must then be rather far above the stellar surface and the situation 
is as discussed by Basko and SunyaevP In contrast, the column solution given by 
Inoue (1975) indicates that the shock would be much closer to the NS surface.^ 
As a limiting case, Braun and Yahel suggested that the shock occurs just on the 
NS surface, and in some cases cannot occur at allP^ However, if the NS boundary 
pressure becomes quite highj^ the shock front may be further out and determining 
its location will depend on complicated physics. A lot of studies of accretion flow 
onto magnetized NSs have considered flow only just near to the NS surface™ global 
treatment of the flow going from the disc to the NS has been considered by only a 
limited number of authorsP*^ 1 ^^ 

The shock location may affect the spectra, luminosities and light curves of X-ray 
pulsars. The result that the shock location may be rather distant from the stellar 
surface has two implications for observations of X-ray pulsars. Firstly, as discussed in 
the previous section, the lower limit for the temperature of the shocked flow (at low 
temperatures for shocked material but much hotter than the unshocked material) 
depends only on the shock location. The region which has reached high temperatures 
due to shock heating, but is still optically thin, corresponds to the transition region 
discussed by Inoue and Basko & SunyaevP^® The spectrum of an X-ray pulsar 
will correspond to the radiation coming from the NS surface, from the optically thick 
radiative region, and from this transition region. The information obtained about 
the transition region is therefore relevant for modelling the observed properties of 
X-ray pulsars. Secondly, if the shock is far above the stellar surface, the radiating 
area presented by the shocked material will be large and the amount of radiation 
coming from the sides of the funnel flow will be significant. This radiation from 
the side-faces ("fan-beams") has a different energy range and pulse phase compared 
with radiation emitted in the direction along the field lines ("pencil-beams"). This 
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difference between the pulse phases makes the observed light curve complicated: 
having a large radiative region due to a distant shock location makes the fan-beam 
component larger, leading to irregular light curves with phase-inversion, twin-peaks, 
peak splitting and so onP'^'DS From observed light curves and spectra, it should 
become possible to deduce the shock locations. 

The analysis carried out in this paper has been very simplified. In actual bright 
X-ray sources, the complicated features which we have neglected may be important, 
especially near to the stellar surface. The main additional features concerned are: 
(i) Radiation pressure: the radiation energy density can become quite high near 
to the stellar surface and radiation pressure can dominate over gas pressure there. 
In Section 16.21 we have estimated the effects of this in just a rough way; we will 
return to give a more thorough treatment of it in future work, (ii) Radiative energy 
losses: in this paper we have neglected energy losses due to radiation emitted from 
the flow. However, this simplification is not valid near to the stellar surface where 
emissivity due to thermal bremsstrahlung becomes significant compared with the 
internal energy density of the flow and then, in the optically thick regime, black 
body radiation will become more efficient. These effects should be included in a more 
complete treatment. The energy losses may make the shock position move nearer to 
the star because they lower the pressure of the column mainly after the shock and 
hence the shock would occur further downstream, (iii) Oblique dipole field: in order 
for a neutron star to be an X-ray pulsar, the dipole axis must not, in fact, coincide 
with the rotation axis, (iv) Relativistic effects: neutron stars are compact objects 
with strong gravitational fields and also the funnel-flow velocity before the shock 
is quite high in some cases. Therefore both general and special relativistic effects 
may change the results quantitatively, (v) Boundary conditions: there are serious 
uncertainties concerning the boundary conditions both at the stellar surface and at 
the point where the funnel flow leaves the accretion disc. All of the above require 
further work but we believe that the present calculations represent a useful first step 
and grasp the overall nature of these accretion flows. 
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